{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import * from utils            \n",
    "\n",
    "tens = [10,20,30,40,50,60,70,80,90,100,110,120,130,140]\n",
    "alpha_all = np.zeros((len(tens), 100))\n",
    "\n",
    "rat_name = 'R'\n",
    "mod_name = '2'\n",
    "sess_name = 'OF'\n",
    "day_name = 'day1'\n",
    "\n",
    "nums = 3\n",
    "angs = np.zeros((len(tens), nums))\n",
    "xlen = np.zeros((len(tens), nums))\n",
    "ylen = np.zeros((len(tens), nums))\n",
    "score1 = np.zeros((len(tens), nums))\n",
    "score2 = np.zeros((len(tens), nums))\n",
    "nn = -1\n",
    "dim = 6\n",
    "ph_classes = [0,1] # Decode the ith most persistent cohomology class\n",
    "num_circ = len(ph_classes)\n",
    "dec_tresh = 0.99\n",
    "metric = 'cosine'\n",
    "maxdim = 1\n",
    "coeff = 47\n",
    "active_times = 15000\n",
    "k = 1000\n",
    "num_times = 5\n",
    "n_points = 1200\n",
    "nbs = 800\n",
    "\n",
    "sspikes_all, xx,yy = get_spikes(rat_name, mod_name, day_name, sess_name, bType = 'pure',\n",
    "                                 bSmooth = True, bSpeed = True)[:3]\n",
    "spikes_all = get_spikes(rat_name, mod_name, day_name, sess_name, bType = 'pure', bSmooth = False, bSpeed = True)[0]\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "\n",
    "num_neurons_all = len(spikes_all[0,:])\n",
    "for nn, n  in enumerate(tens):   \n",
    "    combinations_sample = list(itertools.islice(random_combinations(np.arange(num_neurons_all), n), nums))\n",
    "    for j, comb  in enumerate(combinations_sample):        \n",
    "        sspikes = sspikes_all[:, comb].copy()\n",
    "        spikes = spikes_all[:, comb].copy()\n",
    "        num_neurons = len(sspikes[0,:])            \n",
    "        times_cube = np.arange(0,len(sspikes[:,0]),num_times)\n",
    "        movetimes = np.sort(np.argsort(np.sum(sspikes[times_cube,:],1))[-active_times:])\n",
    "        movetimes = times_cube[movetimes]\n",
    "\n",
    "        dim_red_spikes_move_scaled,__,__ = pca(preprocessing.scale(sspikes[movetimes,:]), dim = dim)\n",
    "        indstemp,dd,fs  = sample_denoising(dim_red_spikes_move_scaled,  k, \n",
    "                                            n_points, 1, metric)\n",
    "        dim_red_spikes_move_scaled = dim_red_spikes_move_scaled[indstemp,:]\n",
    "        X = squareform(pdist(dim_red_spikes_move_scaled, metric))\n",
    "        knn_indices = np.argsort(X)[:, :nbs]\n",
    "        knn_dists = X[np.arange(X.shape[0])[:, None], knn_indices].copy()\n",
    "        sigmas, rhos = smooth_knn_dist(knn_dists, nbs, local_connectivity=0)\n",
    "        rows, cols, vals = compute_membership_strengths(knn_indices, knn_dists, sigmas, rhos)\n",
    "        result = coo_matrix((vals, (rows, cols)), shape=(X.shape[0], X.shape[0]))\n",
    "        result.eliminate_zeros()\n",
    "        transpose = result.transpose()\n",
    "        prod_matrix = result.multiply(transpose)\n",
    "        result = (result + transpose - prod_matrix)\n",
    "        result.eliminate_zeros()\n",
    "        d = result.toarray()\n",
    "        d = -np.log(d)\n",
    "        np.fill_diagonal(d,0)\n",
    "\n",
    "        persistence = ripser(d, maxdim=maxdim, coeff=coeff, do_cocycles= True, distance_matrix = True)    \n",
    "\n",
    "        diagrams = persistence[\"dgms\"] # the multiset describing the lives of the persistence classes\n",
    "        cocycles = persistence[\"cocycles\"][1] # the cocycle representatives for the 1-dim classes\n",
    "        dists_land = persistence[\"dperm2all\"] # the pairwise distance between the points \n",
    "        births1 = diagrams[1][:, 0] #the time of birth for the 1-dim classes\n",
    "        deaths1 = diagrams[1][:, 1] #the time of death for the 1-dim classes\n",
    "        deaths1[np.isinf(deaths1)] = 0\n",
    "        lives1 = deaths1-births1 # the lifetime for the 1-dim classes\n",
    "        iMax = np.argsort(lives1)\n",
    "        coords1 = np.zeros((num_circ, len(indstemp)))\n",
    "        threshold = births1[iMax[-2]] + (deaths1[iMax[-2]] - births1[iMax[-2]])*dec_tresh\n",
    "        for c in ph_classes:\n",
    "            cocycle = cocycles[iMax[-(c+1)]]\n",
    "            coords1[c,:],inds = get_coords(cocycle, threshold, len(indstemp), dists_land, coeff)\n",
    "\n",
    "        num_neurons = len(sspikes[0,:])\n",
    "        centcosall = np.zeros((num_neurons, 2, n_points))\n",
    "        centsinall = np.zeros((num_neurons, 2, n_points))\n",
    "        dspk = preprocessing.scale(sspikes[movetimes[indstemp],:])\n",
    "\n",
    "        for neurid in range(num_neurons):\n",
    "            spktemp = dspk[:, neurid].copy()\n",
    "            centcosall[neurid,:,:] = np.multiply(np.cos(coords1[:, :]*2*np.pi),spktemp)\n",
    "            centsinall[neurid,:,:] = np.multiply(np.sin(coords1[:, :]*2*np.pi),spktemp)\n",
    "        times = np.where(np.sum(spikes>0, 1)>=1)[0]\n",
    "        dspk = preprocessing.scale(sspikes)\n",
    "        sspikes = sspikes[times,:]\n",
    "        dspk = dspk[times,:]\n",
    "\n",
    "        a = np.zeros((len(sspikes[:,0]), 2, num_neurons))\n",
    "        for n in range(num_neurons):\n",
    "            a[:,:,n] = np.multiply(dspk[:,n:n+1],np.sum(centcosall[n,:,:],1))\n",
    "\n",
    "        c = np.zeros((len(sspikes[:,0]), 2, num_neurons))\n",
    "        for n in range(num_neurons):\n",
    "            c[:,:,n] = np.multiply(dspk[:,n:n+1],np.sum(centsinall[n,:,:],1))\n",
    "\n",
    "        mtot2 = np.sum(c,2)\n",
    "        mtot1 = np.sum(a,2)\n",
    "        coords = np.arctan2(mtot2,mtot1)%(2*np.pi)\n",
    "\n",
    "        m1b_1, m2b_1, xedge,yedge = get_ang_hist(coords[:,0], \n",
    "            coords[:,1], xx[times],yy[times])\n",
    "        p1b_1, f1 = fit_sine_wave(m1b_1)\n",
    "        p2b_1, f2 = fit_sine_wave(m2b_1)\n",
    "\n",
    "        x,y = rot_para(p1b_1,p2b_1)\n",
    "        xmin = xedge.min()\n",
    "        xedge -= xmin\n",
    "        xmax = xedge.max()\n",
    "        ymin = yedge.min()\n",
    "        yedge -= ymin\n",
    "        ymax = yedge.max()\n",
    "\n",
    "        angs[nn,j] = ((y[0]-x[0])/(2*np.pi)*360)%360\n",
    "        xlen[nn,j] = 1/x[2]*xmax\n",
    "        ylen[nn,j] = 1/y[2]*ymax\n",
    "        score1[nn,j] = f1\n",
    "        score2[nn,j] = f2\n",
    "        print(angs[nn,j])\n",
    "np.savez('Toroidal_topology_grid_cell_data/Results/toroidal_classification', \n",
    "         angs=angs, xlen=xlen, ylen=ylen, score1=score1, score2=score2)\n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 75,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAOsAAADrCAYAAACICmHVAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjMsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+AADFEAAAH80lEQVR4nO3bva4dVwGG4W9FJMGGKLHkLSgQmBJoT0FJEVHBFdD7Bmi5Bho6XwwSEuXORYQq0g6xE2QnmJ9Fh4Q8ludnrTUz5nnK46P9U7w+a++Zr9RaAxzfO3u/AGAescJJiBVOQqxwEmKFkxArnMS3lvzyw4cP66NHjzq9FOCTTz75vNZ6mfq3RbE+evQo1+u1zasCXlFK+fR1/+YYDCchVjgJscJJiBVOQqxwEou+DQbmu/3t7/n1H/6cj+6/mw/vvZsH99/LDx7cy+9+9dNVjydW6OTZi5f57Ktv8tlX3/z3Zz9++J3VsToGQyfPvv7HKz/78N67qx/vjbGWUh6XUq6llOvtdlv9RPD/5tmLV2N9cL9jrLXWJ7XWu1rr3eUyeRcUMOHpi5ev/Oyj+++tfjzHYOjky4m/rF2PwcA6z76e+ssqVjicp5OfWR2D4XCmjsH+ssIBTR2DfWaFA5q+dOMYDIczFatjMBzQs6nrrPf8ZYVDefnPf+f5y3/9z8/eKckH315/O75YoYPXfbn0zjtl9WOKFTqYvmyz/giciBW6aL24ScQKXbRe3CRihS5aL24SsUIXrRc3iVihi9aLm0Ss0EXrxU0iVuii9eImESt00Xpxk4gVumi9uEnECl20XtwkYoUuWi9uErFCcz0WN4lYobkei5tErNBcj8VNIlZorsfiJhErNNdjcZOIFZrrsbhJxArN9VjcJGKF5qa+Dd5691IiVmhuanGz9e6lRKzQXI/FTSJWaK7H4iYRKzTXY3GTiBWa67G4ScQKzfVY3CRihaZ6LW4SsUJTvRY3iVihqV6Lm0Ss0FSvxU0iVmiq1+ImESs01Wtxk4gVmuq1uEnECk31WtwkYoWmei1uErFCU70WN4lYoalei5tkRqyllMellGsp5Xq73Zo8Kbytei1ukhmx1lqf1Frvaq13l8ulyZPC26rX4iZxDIamei1uErFCMz0XN4lYoZmei5tErNBMz8VNIlZopufiJhErNNNzcZOIFZrpubhJxArN9FzcJGKFZnoubhKxQjM9FzeJWKGZnoubRKzQTM/FTSJWaKbn4iYRKzTTc3GTiBWa6bm4ScQKTfRe3CRihSZ6L24SsUITvRc3iVihianFTcsvlxKxQhOT3wQ3vMaaiBWa6L24ScQKTfRe3CRihSZ6L24SsUITvRc3iVihid6Lm0Ss0MTUMdgXTHBALt3ASfRe3CRihSZ6L24SscJmIxY3iVhhsxGLm0SssNmIxU0iVthsxOImEStsNuKyTSJW2GzE4iYRK2w2YnGTiBU2G7G4ScQKm41Y3CRihc1GLG4SscJmIxY3iVhhM5du4CRGLG4SscJmIxY3iVhhk1GLm0SssMmoxU0iVthk1OImEStsMmpxk4gVNhl12SaZEWsp5XEp5VpKud5uty4vAs5q1OImmRFrrfVJrfWu1np3uVy6vAg4q1G3GiaOwbDJ5K2GHa6xJmKFTUYtbhKxwiaOwXASoxY3iVhhk0NdugFeb9TiJhErbDK5uHEMhmN57eLm/faLm0SssNrIxU0iVlht5OImESusNnJxk4gVVht52SYRK6w2cnGTiBVWG3mrYSJWWG3k4iYRK6w2cnGTiBVWcwyGkxi5uEnECqu5dAMnMXJxk4gVVhu5uEnECquMXtwkYoVVRi9uErHCKqMXN4lYYZXRi5tErLDK6Ms2iVhhldGLm0SssMroWw0TscIqoxc3iVhhldGLm0SssIpjMJzE6MVNIlZY5elzl27gFL50UwScw+jFTSJWWGyPxU0iVlhsj8VNIlZYbOqyzYPOR+BErLDY1OLmw85fLiVihcWePp+61VCscDjTW1bHYDicPW41TMQKi+2xuEnECovtsbhJxAqLOQbDSeyxuEnECovtsbhJZsRaSnlcSrmWUq632637C4Kj22Nxk8yItdb6pNZ6V2u9u1wu3V8QHN0ei5vEMRgW2Wtxk4gVFtlrcZOIFRbZa3GTiBUW2Wtxk4gVFtlrcZOIFRbZa3GTiBUW2etWw0SssMhei5tErLDIXoubRKywiGMwnMRei5tErLDIXoubRKywyF6Lm0SssMhei5tErDDbnoubRKww256Lm0SsMNuei5tErDDbnoubRKww256Lm0SsMNvUX1bHYDigqc+sjsFwQHsubhKxwmx7Lm4SscJsey5uErHCbHsubhKxwmx7Lm4SscJsU4sbl27ggKYWNy7dwMHsvbhJxAqz7L24ScQKs+y9uEnECrPsvbhJxAqz7L24ScQKs+y9uEnECrPsvbhJxAqz7L24ScQKs+y9uEnECrPsvbhJxAqz7L24ScQKs+y9uEnECrPsvbhJxAqz7L24ScQKb3SExU0iVnijIyxukqTZfw1fPH+ZWmurh4PD+MtfX7zys9GfV5OGsX78+z/li4mbneFtNPrzauIYDKuMvmyTzIi1lPK4lHItpVxvt9uI1wSHd/ng/eHP+cZYa61Paq13tda7y+Uy4jXB4X38k+8Nf85mn1lH3ycJe7h89/385uc/zC9/9v3hz90s1j/+9hetHgqY4AsmOAmxwkmIFU5CrHASYoWTECucRFly830p5Zbk0xm/+jDJ52tf1EF4D8fxNryPue/hR7XWybuPFsU6VynlWmu9a/7AA3kPx/E2vI8W78ExGE5CrHASvWJ90ulxR/IejuNteB+b30OXz6xAe47BcBJihZMQK5yEWOEkxAon8R/TvxltHVaWfAAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "tens = [10, 20, 30, 40, 50, 60, 70, 80, 90,100,110,120,130,140]\n",
    "\n",
    "f = np.load('Toroidal_topology_grid_cell_data/Results/toroidal_classification.npz', allow_pickle = True)\n",
    "angs = f['angs']\n",
    "xlen = f['xlen']\n",
    "ylen = f['ylen']\n",
    "score1 = f['score1']\n",
    "score2 = f['score2']\n",
    "f.close()\n",
    "\n",
    "plt.figure()\n",
    "ax = plt.axes()\n",
    "successes = np.zeros(len(tens))\n",
    "fit_thresh = 0.25\n",
    "for i in range(len(tens)):\n",
    "    successes[i] = np.sum(((angs[i,:]<70) & \n",
    "                           (angs[i,:]>50) &\n",
    "                           (np.abs(1-xlen[i,:]/ylen[i,:])<0.25) &\n",
    "                           (score1[i,:]<fit_thresh) & \n",
    "                           (score2[i,:]<fit_thresh)))\n",
    "    \n",
    "ax.plot(np.concatenate(([0], np.divide(tens,140)))*100,\n",
    "        np.concatenate(([0], successes/nums*100)), label='0.25', lw = 5)\n",
    "ax.set_xticks(np.divide([0,35,70,105,140],140)*100)\n",
    "ax.set_xticklabels(('','','','',''))\n",
    "ax.set_yticks((0,25,50,75,100))\n",
    "ax.set_yticklabels(('','','','',''))\n",
    "ax.set_aspect('equal', 'box')"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}
